Semi-classical approximation for the second harmonic generation in nanoparticles 
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Second harmonic generation by spherical nanoparticles is a non-local optical process that can also be viewed 
as the result of the non-linear response of the thin interface layer. The classical electrodynamic description, 
based e.g. on the non-linear Mic theory, entails the knowledge of the dielectric function and the surface non- 
linear optical susceptibility, both quantities are usually assumed to be predetermined, for instance from experi- 
ment. We propose here an approach based on the semi-classical approximation for the quantum sum-over-states 
expression that allows to capture the second-order optical process from first principles. A key input is the elec- 
tronic density, which can be obtained from effective single particle approaches such as the density-functional 
theory in the local density implementation. We show that the resulting integral equations can be solved very 
efficiently rendering thus the treatment of macroscopic systems. As an illustration we present numerical results 
for the magic Na^g^^ cluster. 
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I. INTRODUCTION 

For the description of a wide range of phenomena in physics 
and chemistry one is faced with the question of how to predict 
in a reliable and system specific way the response to external 
electromagnetic fields that impatt energy co and possibly mo- 
mentum k to the system [1]. To name but few examples, we 
mention here the response of nanoparticles to light in the lin- 
ear [2] and the non-linear regimes [3]. For electrons as a per- 
turbation we refer to the overview [4]. A well-studied route 
to address these issues is the concept of the dielectric response 
functions which is detailed in standard books such as in [1], 
but also in recent reviews, e.g. in Refs. [2-4]. 

For extended bulk metals and metalUc surfaces the dielec- 
tric response was extensively studied (Ref. [5, 6] and refer- 
ences therein). An extension to finite-size systems brings 
about a number of new aspects that complicate the theory 
but at the same time offer new opportunities for the occur- 
rence of new phenomena. For example, it is known that 
the second harmonic generation (SHG) is forbidden in sys- 
tems with an inversion symmetry on the dipole approxima- 
tion level. Therefore Ostling, Stampfli, and Bennemann [7] 
and Dewitz, Hiibner, and Bennemann [8] proposed the an- 
harmonic oscillator model to describe the second-order non- 
linear response from spherical systems. A further step was un- 
dertaken by Dadap et al. [9] who considered the small-particle 
limit {ka < 1) and described SHG as a mixture of dipole and 
quadrupole excitation processes {a stands for the system size). 
It was also shown how to connect the model output to exper- 
imentally observed quantities, i.e. the elements of the non- 
linear optical susceptibility tensor )^^\oj). In this approach 
the inversion symmetry is not a necessary assumption as it is 
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assumed that the second harmonic generation originates from 
the thin surface layer, where the symmetry is broken anyway. 
This idea can be extended in several directions as was demon- 
strated in numerous works: Pavlyukh and Hiibner [10] devel- 
oped the non-linear Mie theory valid for particles of arbitrary 
sizes, de Beer and Roke included the sum-frequency gener- 
ation mechanism into the considerations [11], the cylindrical 
geometry was treated by Dadap [12] and finally the theory for 
arbitrarily shaped particles was developed by de Beer, Roke, 
and Dadap [13]. 

All these theories based on classical electrodynamics rely 
on the knowledge of the frequency-dependent dielectric func- 
tion and the non-linear optical susceptibility tensor. With the 
fabrication processes being perfected and the system sizes 
lending smaller and smaller one may wonder to which ex- 
tent quantum effects are important and whether it is justified 
to use the same susceptibility tensor to describe semi-infinite 
and finite size systems on the nanometer range. To shed light 
on these issues, it is desirable to have a quantimi theory for 
the non-linear response on the nanoscale. The fully atomistic 
approach seems to be out of reach for present computers, as 
currently maximally hundreds of atoms are possible to treat 
using quantum chemistry codes. Yet, the outstanding question 
is, how important are the electronic correlation effects and is 
there possibly a way to stay on the solid quantum theory basis 
while treating larger systems? 

There is an affirmative answer to these questions as was 
demonstrated recently in the linear optics case by Prodan and 
Nordlander [14]. They succeeded to push the limits of the 
time-dependent density functional theory (TDDFT) to metal- 
lic systems containing millions of atoms. But at the same 
time they demonstrated that for these system sizes the semi- 
classical approach becomes very accurate. This is a marked 
finding as it allows to replace the complicated sum-over-states 
quantum mechanical expression [15] for the polarizability ten- 
sor with a single integral equation. Consequently, there is 
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only one parameter in the theory: the ionic density distri- 
bution. With some reasonable assumptions about the ionic 
density such as in the jellium model (this assumption is rea- 
sonable even for molecular structures, as we have shown re- 
cently [16, 17] for fullerenes. The usefulness of the jelUum 
model was demonstrated by the pioneering works of Ekardt 
on sodium clusters [18, 19] or of Puska and Nieminen on 
Cfio [20].) we can obtain the ground state electronic density 
from the solution of the Kohn-Sham equations and express the 
response function in its terms. The Drude dielectric function 
follows automatically. 

The semi-classical approximation is rooted in the work 
of Mukhopadhyay and Lundqvist [21] who derived the cor- 
responding integral equation within linear response theory. 
Their theory was applied in numerous cases ranging from 
plasmons in the jellium model to collective resonances in 
Cgo [22] or carbon nanotubes [23]. The equation can also 
be derived starting from the quantum mechanical sum-over- 
states expressions [24] and using the assumption that the fre- 
quency of the external field is large compared with the single- 
particle gap {oj » Eg) [25]. 

One may wonder why this program was not implemented 
for nonlinear optical processes. As a matter of the fact, al- 
ready in 1973 Wang, Chen and Bower [26] classically treated 
second harmonic generation from alkali metals. A decade 
later Apell [27] derived an expression for the second harmonic 
current in the form: 

j(r, 2co) = a[/(E • V)E + |3(V/)£2 + y/(V • E)E], (1) 

where for the unperturbed ground state electron density in the 
form n{r) = no fir) the a parameter is a function of no and 
o), whereas |3 and y are constants. The theory gained less at- 
tention than, for example, the Mukhopadhyay and Lundqvist 
work because the connection to quantum mechanics was lost. 
Here we mention nonetheless numerous works in the field 
extending over more than three decades: the small homoge- 
neous spherical particle limit [28], the Rayleigh-Gans scatter- 
ing approximation for a lattice of such particles [29], second 
harmonic generation by two-dimensional particles [30], or a 
more recent treatment of arbitrary geometries [31]. Notice 
that the assumption of homogeneous electron density distri- 
bution within the sample is an inevitable component of such 
classical theories. 

It is possible to revive the theory by noting that the electric 
field (E) in Eq. (1) must also include the induced field. This 
simple observation immediately raises the level of the theory 
to the random phase approximation (RPA) or, if we include 
electronic correlations, to TDDFT level. Our manuscript is, 
hence, a formalization of this message. 

To this end we first derive an expression analogous to (1) 
starting from the sum-over-state quantum mechanical formula 
for the non-linear optical susceptibihty [15] and employing 
the high-frequency approximation. Although quite technical, 
we believe that this derivation (Appendix A) has its own mer- 
its as it estabhshes the equivalence of the hydrodynamic ap- 
proximation of Apell and the high-frequency semi-classical 
expansion. It is interesting to notice that the asymptotic 
behavior of this quantity is at variance with the results of Scan- 



dolo and Bassani [32] who predict a l/oj^ decay. This seems, 
however, to be a direct consequence of the assumption of the 
aU-dipole excitation mechanism that underly their study. It is 
possible to prove from general principles that the inclusion of 
the quadrupole excitations leads to the l/co'^ asymptotics [33]. 

Our derivation raises the question of whether it is sufficient 
to know the unperturbed ground state density to obtain the 
lowest order approximation for an arbitrary response function. 
We recall that from the point of view of the diagrammatic per- 
turbation theory [34] SHG comprises three processes in which 
the 2a) photon is emitted before, between or after the absorp- 
tion of two w-photons. Consequently, one might wonder if 
each diagram of this expansion can also be expressed in terms 
of n(r). As we demonstrate below, the answer is negative, 
one additionally needs the one-particle density matrix y(r, r') 
whose diagonal elements are given by n(r). This comes not 
as a surprise if we consider the analogy with the orbital-free 
kinetic density functional theory [35] where this matrix enters 
the kinetic energy term. 

The non-linear susceptibility relates the second-order in- 
duced density to the local electric field. The latter is a mi- 
croscopic quantity that can be connected to the external field 
by knowing the linear response function. In the linear case 
it is a standard route to get the RPA dielectric response. In 
the non-Unear case the procedure is, probably, less known. 
Therefore, we follow here a very pedagogical treatment of 
Liebsch and Schaich [36]. In fact, they apphed a trick sug- 
gested by Zangwill and Soven [37] to avoid the sunmiation 
over the infinite number of unoccupied states for the calcula- 
tion of the response functions. This allowed them to describe 
SHG at simple metal surfaces as effectively one dimensional 
systems (it is basically the same approach that enabled Prodan 
and Nordlander [14] to treat very large spherical systems). 

The outline of this work is as follows: In Sec. II we start 
with the most general relation between the induced densi- 
ties and the effective fields [Eq. (2.5) of Liebsch and Schaich, 
Phys. Rev. B 40, 5401 (1989)] and formulate integral equa- 
tions for the induced density with the non-interacting response 
functions as kernels. In Sec. Ill we discuss the case of a spher- 
ical symmetry and the simplifications it implies for the numer- 
ics. Finally, the second harmonic response of the magic Najggg 
cluster is presented in Sec. IV for the illustration of our the- 
ory. Based on our recently developed method for the solution 
of integral equations of this type we are able to drastically re- 
duce the computational cost from 0(N^) to 0(N), where A'^ is 
the number of mesh points to represent the density. 

We use atomic units, i. e., ?! = e = m^. = 4nso - 1 through- 
out. Two appendices contain mathematical details to make the 
exposition in Sees. II-IV self-contained. 



II. INTEGRAL EQUATIONS 

Our theory can easily be extended to include electron cor- 
relation effects by using the exchange correlation functional 
of DPT. Our formulation here is at the random phase approx- 
imation level. Within this approximation the density-density 
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response function can be obtained as 

xu{ri -r2)x(r2,r'; CO), (2) 

where i;(r-ri) is the Coulomb potential and;^f^"'(r, r';oj) is the 
non-interacting density-density response function (known as 
Lindhard function for the homogeneous electron gas in three 
dimensions, 3D): 



y fi~ fj 

l^uj + Ei-Ej + irj 

X (^,(r)^}(r)^^/r')(^;(r'), 



where / is the Fermi function and j refer to the collections 
of quantum numbers that uniquely characterize the electronic 
states of the system. The infinitesimally small positive num- 
ber rj shifts the poles from the real axis and ensures, thus, the 
causaUty of the response function. In what follows we wiU 
assume it can be incorporated in the oj variable. 

Let us consider the response of the system subject to the 
harmonic electric field oscillating at the frequency co, i.e. 
^''''(r; t) = (p''^\r) e^"^'. Then, the induced density which os- 
cillates at the frequency of the appUed field is given by: 

<5n(')(r) = J dr'x{r,r'-u)ip\r') 

= r<irVO)(r,r';a>)^(i\r'), (4) 



where ^"^''(r) is the induced local field oscillating at the funda- 
mental frequency and consisting of the external potential plus 
the Hartree potential corresponding to the induced density: 

<p^^\v) = ip^^\Y) + ^dr'v{v-r')6n^'^\Y'). (5) 

The induced density oscillating at the double frequency results 
from the non-linear process described by the X2^ response 
function and from the Unear response to the local field (p^^^(x) 
oscillating at Ico: 

6n^^\v) = j ir'j ir"xf\r; r', r"; w) (p^^\r') (p'-^\r") 

+ f dr'x'-'^\r,r';co)ip^^\r'). (6) 



Because there is no external field at 2a> the local field at this 
frequency is given by the Hartree potential: 



while eqs. (6) and (7) result in the integral equation for the 
second harmonic density: 

6n^^\r) = f^^>(r)+ 

J ir'J ir"x^'^\r, r'; w) y(r' - r") 6n^^\r"). (9) 

We defined the source terms following [36] as: 

^(i)(r) = Jdr'/"''(r,r';Lo)(p^"\r'), (10) 
^(2)(r) = J ir'J ir"xf\r; r', r"; co) (p^^\r') <p''^\r"). (1 1) 

Eqs. (8, 10) and eqs. (9, 11) when coupled with an appro- 
priate approximation for the non-interacting response func- 
tions allow to completely describe the linear and the second- 
harmonic response. 

The semi-classical approximation for the ^^\r) has already 
been derived previously [24] with the result: 

^(i)(r) = ^(Vn(r) • V^(°'(r) - n(r)A^(°'(r)). (12) 

In Appendix A we derive the second-harmonic generation 
source term (A20) that can be represented as: 

^''(r) = ^ V ■ |v^(i'(r) [n(r)A^(»(r) + (v^(i>(r) • Vn(r))l 

+ ^n(r)V(v^(i>(r))']. (13) 



Although these expressions look rather complicated they can 
be further simplified in the case of spherical symmetry. 



III. SPHERICAL SYMMETRY 

For spherically symmetric systems (i. e. n(r) = n(r)) 
eqs. (8) and (9) simphfy considerably when projected onto 
spherical harmonics. Since these equations have the same 
functional form we introduce an index / - 0, 1, 2 to label 
corresponding quantities. We use the multipole expansion of 
the fields: 



(m 

This is a general form consistent with the spherical symme- 
try. Similar expansions will be used for the densities and the 
source terms: 



Equations (4) and (5) yield the integral equation for the linear 
density: 

6n'-\r) = f*"(r)+ 

J ir'J ir"x^''\r, r'; (o) y(r' - r") 6n^^\r"), (8) 



For numerical calculations in this work we restrict ourselves 
to the mixture of the dipole and quadrupole excitations (ipf^ = 
for ^ > 2). Further simplifications can be achieved by con- 
sidering the external field to be a plane- wave: 

exp(jkr) =4nY, i' Mkr)Y(„{Sl,)Yf„{a,), (14) 

tm 
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and when we align the coordinate system along the k- 
direction: 



exp(*r) = = AnY,i' ^^^Mkr)Yeo (^) • (15) 



When the wave-number k is small compared to the dimen- 
sion of the system a (i. e. fea < 1) it is justified to replace 

the spherical bessel functions with their small argument ap- 
proximations jtiz) ~ z" l(2{ + 1)!!, use the definition of the 
multipole moments: 



and express the external potential as: 

^^'^(r) = 27#4TG»(rX or 



(2^-1)!! 



(0)/ X iik) 



(2^-1)!! V 2^ + 1 



An 



(16) 



(17) 



where each term is a harmonic function i. e. A(2fm(r) = 0. 
Thus, in this approximation, it is sufficient to consider only 
the first term in (12): 



f(i)(r) = — Vn(r)- V^^"'(r). 



1 



,(0)/ 



0)^ 



(18) 



With the help of eqs. (4) and (5) we have: 

nA^'i' = -0}], 6n^^\ (V^'^' • Vn) = {o? - w^) 6n^^\ 

where we introduced the local plasmon frequency in analogy 
with the expression for homogenous systems ^^(r) = Ann{r). 
Finally, the second-order source term can be given as a sum 
of two contributions ^^\y) = ^^\r) + ^^'>\r): 

= • (6n('\r)[o/ - 2w2(r)) V^("(r)j,(19a) 

^'*'(r) = ^ V • |n(r)V (v^('>(r))' j. (19b) 

The physical meaning of these two terms can be inferred 
by introducing the induced electric field E - -V(^*'' in e.g. 
Eq. (13). The first term contains a contribution from the elec- 
tric quadrupole term E(V • E) and from the density variation 
(the surface dipole term) E(E ■ V«), whereas the second term 
upon the use of the vector identity j = Ex VxE-i-(E- V)E 
and one of the Maxwell equations E x V X E = ioi/cE x H can 
be interpreted as the magnetic "bulk" term and a surface term. 
According to the analysis of Wang, Cheng and Bower [26] and 
ApeU [27] the contribution from f '^''^ is small in the case of the 
non-linear response from the surfaces and the incident elec- 
tric field polarized in the plane of incidence. These arguments 
loose their power in the case of the non-linear response from 
spherical objects. Thus, both terms must be considered. We 
will show below that the treatment of the second term is much 
more involved. Therefore, to illustrate our approach only the 
first term will be included into the numerical algorithm. 



A. Computational scheme 

With these results in hands the numerical algorithm can be 
formulated as follows: 

i) The integral equation (8) is solved with a source term 
(18). It yields the first order density 6n^^^(r) and, there- 
fore, the local potential <p^^^{v); 

ii) Eq. (19) is applied to generate the source term for the 
equation (9); 

iii) This integral equation is solved similarly to (8) in order 
to obtain the second order density. 



B. Linear response 

In this subsection we focus on the first point of the program. 
Already on this level our approach is valuable as it provides 
the optical absorption cross-section. 

a. The linear source term: Using (18) and the explicit 
form of the potential (17) we obtain: 



4>) = .n'(04^, 



(20) 



where n'(r) denotes the derivative of the equilibrium density 
function with respect to the radial coordinate. 

b. The integral kernel: The response functions are in- 
variant under the rotations of the system as a whole. Thus, in 
the linear case it can be expanded as: 

Ar(r,r';w) = YjXi,{r,/;oS)Yl{mi,m. 
In 

Thus, the integral kernel of eqs. (8) and (9) can be projected 
onto the angular momentum eigenstates: 

Kt{r,/) = ^d£lYlja) 

X JdQ'V«(r)-Vr(;(|r-r'|)7^„(Q'). (21) 

It can be integrated using the spherical harmonics expression 
of the Coulomb potential [see Sec. 3.6 of Jackson (Ref. 38)]: 

^2X2^^1 |r>';(")W^ 

where r< (r>) is the smaller (larger) of r and r'. The gradient 
of the spherical harmonics need not to be considered in view 
of the symmetry of the density, the angular integration can be 
done beforehand and we obtain for (21): 



Kt{r,r')^ 



27TT"^'"^^i;frj 



An 



2i+\ "'r'M^t 



n(r)^Gf\r,r'), (22) 
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with the spherical ^-pole Green function defined as 

Gf^(r,r') = ({+l)e(r-r')i'-\ - {0(r' - r). (23) 

With the help of (22) the integral equations (8) and (9) can be 
written in the unified form: 



1 - 



^finio)^ 47r n'(r) 
2^+l~^' 



where the density is obtained from the solution of (24). It is 
easy to see that in accordance with the earlier assumption it 
is sufficient to restrict ourselves to the cases of i = 1,2 and 
m = 0. The most difficult part of the derivation: the transi- 
tion from v?*''(r) to ^*^'(r) can probably be obtained in closed 
form, however, for practical applications it is sufficient to have 
shorter form small-^ solutions. They can be obtained with the 
MAPLE computer algebra package: 



J d>' Gl {r,r)6nYJr ■,lS). (24) 2V5;rL ^ 



This equation is central for our theory in both the linear and 
the non-linear cases. An efficient method of its solution ex- 
ists [24]. In the linear case the equation takes a more sym- 
metric form when re-formulated in terms of the polarizabUity 
function. 

c. Observables: The. C-pole frequency -dependent polar- 
izability deffned as: 

atmiio) = - Jdr J ir' Qf^ir)xir, r'; w)2^m(r'), 

can be computed as the response to the held (f>f\r) — 
— Let us introduce the position-dependent polarizabil- 




adr; aS) = 



4n 
2(+l 



6n^l>ir,oj) 



and use the explicit form (20) for the linear source. Substitut- 
ing these deffnitions in (24) we obtain the following integral 
equation: 



a/(r; a>) - af\r, a>) 



(25) 



where 



(0)/ \ 

a) (r; w) 



47r r"'n'(r) 



In the case of the dipolar response our result (24) coincides 
with Eq. (5) of Prodan and Nordlander (Ref. 14). Finally the 
frequency-dependent polarizabiUty is represented as the inte- 
gral: 

r*oo 

am{o))= I dr/^^af\r;oS). (26) 
Jo 



C. The non-linear source term 

The source term for the second-order response is consider- 
ably more complicated. The central quantity is the induced 
potential ^^'^(r). It can be evaluated by the integration of (5): 



,(0),- 



2^+1 J ri- 



,,.i<5«S('-'). (27) 



?2 



1 



2V57r 



-I- u {n'if{ + n'2f[ - STrnini) 



f[(u'ni + Mn'i) - -j[nifi + 4nr^n{j 



, (28) 



(29) 



where we introduced for brevity = ^^^\r), u = a/ - 
2colir), nt = 6n^^^{r\ ft = <p^^{r). 

The second part f ^^'''(r) is presented in Appendix B for ref- 
erence, however, will not be included into the numerical algo- 
rithm. 



IV. NUMERICAL RESULTS 

In the present contribution we continue to study the optical 
properties of the magic Najggg cluster (Fig. la). This sys- 
tem simultaneously possesses completely closed geometric 
and electronic shells. This makes it exceptionally stable and 
attractive for the numerical calculations: its electronic struc- 
ture can be obtained easily by using the density functional ap- 
proach. We do not pursue here a fully atomistic approach, 
but rather solve the radial Kohn-Sham equation (using the 
renormahzed Numerov method) in the presence of the spheri- 
cally symmetric ionic density. We either use the standard jel- 
lium model which starts from the homogeneous ionic density 
(Fig. lb) or the density is obtained from the realistic geometric 
model by applying the gaussian broadening of ~ 1 A width to 
ionic positions and performing a spherical averaging (Fig. Ic). 
In accordance with the proposed numerical scheme we present 
here a) the dipolar and the quadrupole Unear optical absorp- 
tion coefficients afX){(jL)) and the induced densities 5nj,Q(r;w) 
{( — 1,2) (Fig. 2); b) results for the second-harmonic source 
^'•^'(r); and c) the non-linear dipolar and quadrupole optical 
responses at 2oj frequency (Fig. 3). 

While the linear optical response of metalUc clusters is 
well understood: typically the optical absorption proffle only 
sUghtly deviates from that of the idealistic system: a sphere 
filled with an electron gas of constant density. There, the opti- 
cal absorption coefficient is peaked at the energies of the sur- 
face plasmon modes: 



2€+\ 
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10 20 30 

Radial distance ;■ (A) 

FIG. 1. (Color online) a) Geometrical structure of 10-shell icosahe- 
dral Na^ggj cluster. The next two panels show the electronic structure 
(only £ = states are shown), the Kohn-Sham potential (black solid 
line), the electron density (red solid line) from the self-consistent lo- 
cal density approximation (LDA) calculations for the jellium model, 
the dotted line marks the ideal jellium background (b) and for the 
spherically averaged realistic ionic density (dashed blue line) (c). 
The intershell spacing is fixed at the bulk nearest neighbor distance. 



The position of the maxima is only weakly dependent on the 
details of the electronic structure: our calculations almost per- 
fectly match corresponding idealistic values cJi - 3.45 eV and 
t02 = 3.78 eV for the bulk Na density (r, = 3.96). The spill- 
off of the electron density in realistic systems mostly leads to 
the broadening of the surface plasmon resonances. To illus- 
trate this fact we choose the off-resonance value of the fre- 

(2) 

quency (w* — 5.75 eV) and plot the source ^J, (r, a>*) and the 

(2) 

induced density dn^ (r, oj*) (Fig. 2). While in the idealistic 
situation the induced density is distributed over the surface 
of the sphere (where the electron density abruptly changes) 
in the realistic case we observe numerous features associ- 
ated with slow electronic density variations within the cluster. 
They contribute to the optical absorption in the off-resonance 
regime. However, the relative weight of these oscillations de- 
creases when the frequency approaches the resonance. There, 
the fast density variation at the surface dominates the spec- 



trum. 

The non-linear optical properties (Fig. 3) of even such 
simple systems are not fully understood. It is interest- 
ing to observe that two different excitation mechanisms (the 
quadrupole transition at co or 2a> frequency) lead to almost 
identical frequency dependence. Unlike in the linear case, the 
efficiency of the frequency conversion vanishes at the plas- 
mon resonance and has two pronounced peaks at the energy 
slightly below and above. We also do not observe a strong 
correlation between the source and the induced density as in 
the off-resonant linear case (cf. blue and green curves). How- 
ever, the spatial variation of these quantities is not erratic (as 
the plots might be suggesting). The complicated radial de- 
pendence is the result of the derivatives of the induced den- 
sity and the potential in the non-linear source terms. Thanks 
to the linear-scaling algorithms we are using at each stage of 
the computation we are able to eliminate any spurious contri- 
butions from the numerical differentiation while ensuring the 
convergence of our results with respect to the number of dis- 
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FIG. 2. (Color online) The linear ( = I (top row) and f = 2 (bot- 
tom row) optical response of the icosahedral Na^g^^ cluster based on 
the standard jellium model (black solid line) and on the model with 
the realistic ionic density (red dashed line). Right panels show the 
spatial dependence of the corresponding source terms ^f \r,CL)') and 
the resulting induced densities Sn f\r, o)') for a particular frequency 
value oj' = 5.75 eV. 
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FIG. 3. (Color online) The second-order non-linear { = I (top row) 
and £ = 2 (bottom row) optical response of the icosahedral Najg^^ 
cluster for the two models at Fig. 2. Right panels show the spatial 
dependence of the corresponding source terms ^^\r, w*) and the re- 
sulting induced densities 6nf\r, a)') for a particular frequency value 
w* = 5.75 eV. 



cretization points and the value of the broadening parameter 
'7- 



ments allowing, thus, for the treatment of larger systems. The 
work along these lines is already in progress. 

We believe that our method is sufficiently versatile as to 
allow for further extensions. In particular, it is straightfor- 
ward to modify the method for systems with axial symmetry, 
or even to consider the symmetry-free case. To treat larger 
systems one must include higher multipole moments. This 
poses a question of how to find the second-order source term 
in this case. We believe that a direct manipulation with the 
MAPLE computer algebra system rather than a formal solution 
in terms of 3 j symbols is feasible. 

Regarding the physical aspects of our approach: The fact 
that it is free of any adjustable parameters is not necessar- 
ily beneficial. In fact, for metallic systems with localized 
(i-electrons peculiarities of the electronic structure might be 
reflected in the optical properties. In this case the classical 
electrodynamics approach with the experimentally measured 
dielectric function and susceptibility tensor might give more 
accurate results. On the other hand, for systems with simpler 
electronic stiTicture our method is capable of taking into ac- 
count quantum effects such as the spill-out of the electronic 
density. 

Finally, our work establishes the equivalence between the 
semi-classical approximation and the hydrodynamic approach 
of Apell [27]. The latter, however, is just a classical theory if 
not corrected for local field effects. We also touched upon the 
question of the representability of the response functions in 
terms of electronic densities and show that in general a more 
complicated quantity such as the one-particle density matrix 
must be introduced. However, for the second-harmonic gen- 
eration these terms cancel in the final expression. 



ACKNOWLEDGMENTS 

The work is supported by DFG-SFB762 and DFG-SFB/TR 
88. W. H. gratefully acknowledges the hospitality of the 
"Nonequilibrium Many-Body Systems" group at Martin- 
Luther-University Halle- Wittenberg during his sabbatical. 



V. CONCLUSIONS 



Appendix A: Semi-classical expression for the first 
hyperpolarizability 



We developed a semi-classical theory of the second- 
harmonic generation in spherical particles. Although it origi- 
nates from the exact quantum-mechanical sum-over-states ex- 
pression and takes the local fields into account it is free from 
the small system size limitation. Because an efficient method 
for the solution of the corresponding integral equations ex- 
ists the theory can be applied even to macroscopic systems 
provided the electronic density can be found. For this pur- 
pose the jellium approximation for the ionic density can be 
used as we illustrate by computing the SHG spectrum of the 
Najggci cluster The scattering cross-section and the angular 
distribution of the intensity are other important experimental 
observables. Calculation of these quantities will be presented 
elsewhere together with the inclusion of higher multipole mo- 



We start with the expression for the second order non-linear 
response: 

<A:(r).A,(r)^*(r')«A/r')-A*(r")«A*(r") 
4 '(r; r ,r ; w) = > ^ ^— "t^ ^r^. 



ij,k 
fk - fj 



cl> + Ej - Eit + irj ci) + Ei - Ej + irj j 



(Al) 



where / is the Fermi function and /, j, k refer to collections of 
quantum numbers that uniquely characterize electronic states 
of the system. In what follows we use the the following no- 
tations Ei - E/i - Ell-, fi — fk - fik, etc. The infinitesimally 
small positive number 77 shifts the poles from the real axis and 
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assures, thus, the causality of the response function. In what 
follows we will assume that it can be incorporated in the co 
variable. 

Let us consider a generic type of integrals: 

<D(r) = J d{r'Y")xf{r; r', r"; cS) 0(r') 0(r") (A2) 

and introduce our basic approximation. The semi-classical 
approximation (SCA) implies a high frequency condition \Ei— 
Ej\ «: cj. Thus, the above expression can be expanded as 
a power series of o). The term proportional to l/w^ in the 
expansion of <l)(r) vanishes in view of the completeness of the 
electron wave- functions: 



2^/';(riM,(r2) = 5(ri-r2). 



(A3) 



In what follows we will assume all wave-functions to be real. 
This can always be done without the loss of generality for 
time-invariant systems. Odd power terms vanish in view of 
the symmetry of the expression with respect to permutation 
i k and r' <z> r". 

The 1 Ilo^ term has the following form: 

<b^\v) = f dir'r") 2 Xij,,(r, r', r") 0(r') 0(r") 



1 



(A4) 



The first term vanishes after using the property (A3), the in- 
tegration over r', and exploiting the permutation symmetry 
j ^ k and r <:± r" of the expression under the integral. The 
term proportional to can be re-written as lEn^Ej^fj and 
combined with the second terms. Finally our expression can 
be written as: 

a)('*)(r) = io^^-^'Cr) - ^^*''\r) + ^^^*'\r) (A5) 

with following notations: 

¥^"\r) = f dir'r") ^ X,-,,-,,(r, r', r")[Elfk\4>(r') 0(r"), 

^ i,j,k 

¥^''\r) = f dir'r") J] X^j^^ir, r', r")[£»£;,/,l0(r') 0(r"), 

^ i,j,k 

<D(4^)(r) = f dir'r") ^ Xtj,kir, r', r")[£,vt£;i/*]0(r') 0(r"). 



Uj,k 



In what follows we will make use of 

2 



Eij^iir)il,jir)^--V -^ijir). 



(A6) 



where ^y(r) - i//jir)Vi//iir) - i/fi(r)V^/r). This foUows from 
the Schrodinger equation. 

Let us consider first the (S)^'^''\r) term. After the summation 
over j and integration over r" we arrive at: 

^^'"Kr) = \ f dr'<l>\r') ^ /,[V, • ^,,(r)][V^ • ^,,(r')]. 

^ i,k 



We further pull out of the integration and use the Gauss 
theorem to apply V,- to the (jP'ir'): 

a,(4a) = _Z: . J dr' 2 A^*(r) [^,-,(r') • Vr><t>\r')l (A7) 



Now the summation over i and k can be performed by intro- 
ducing the one-particle density matrix: 



^fiil/iiri)tlfiir2) = r(ri,r2), 



(A8) 



with the diagonal elements given by the density: 

2/;-^Kr)^,(r) = n(r), (A9) 

i 

Quite general we can represent the density matrix in the form: 

r(ri,r2) = /(ri)/(r2)6((|ri - r2|). 

From the definition (AS) a useful integration formula can be 
derived: 

J" dr'6ir - r')F(r')V,y(r, r') 

= ^ dr'6ir-r')Fir')Vr'yir,r') = ^Fir)V Mr). (AlO) 

Let us now retum to (A7) and perform the sum over the states 
using eqs. (AS) and (A9). 

^^^"^ = Z r dr'Vr"pHr') 

ab ^ 

X [7(r,r')V;.'V|!,(5(r,r') - V;,'r(r, r')Vf,5(r, r') 
- V^5(r, r')V* r(r, r') + 5(r, r')V^Vf,7(r, r')]. (All) 

In order to evaluate the integral we must remove all diff'eren- 
tial operators from the (5-function. 

O^^"' = - ^ 2] K r dr'V^^(p\r') 

ab ^ 

X [V* V;:'[(5(r,r')r(r,r')] - 2V«[(5(r,r')V*r(r,r')] 
- 2Vf, [6iY, r')V°y(r, r')] + 45(r, r')V^Vf,7(r, r')]. (A12) 

It can be further simplified by using Green's first identity for 
the first and third terms: 

^(4a) = ^ J" dr'5(r,r')[A(^^(r')r(r.r') 
-h2V*02(r')V*7(r,r')l 
" ^Z^" / ^r'^(r>r')[A[0'(r')]V^r(r,r') 

+ 22]v«V*y(r,r')V^02(r')l. (A13) 



Let us introduce the tensor T'''{r) = limr'_,r V^V^-y(r, r') and The 0^'**''^(r) term requires special attention: 
integrate: 



- i V ■ [Vn Acf,^] - V''[r''vV2j. (A14) 

ab 



Using (A6) we re-write O^"'*' in the form without explicit 
energies: 



<l)<4'')(r) = ^ 



i,j,k 

X iffiir')iffjir'). (A15) 



We perform the summations over the states and use the Gauss 
theorem to transfer V^" on ^(r") to obtain 



3,(4fe)(r) = - J 2 V« r J(r'r")0(r') V*,0(r") 

ah ^ 

X [6(r, r")V^5(r, r') V';„r(r", r') 

- Vf„(J(r, r") V^5(r, r') r(r", r') 

- V^(5(r, r") 5(r, r') Vf„y(r", r') 

+ Vf„ V^(5(r, r")) 5(r, r') r(r", r')] (A16) 



The integration yields: 



= -\YjK[<P(rW''r [ ^/r"V^,0(r")V^,[5(r,r")r(r",r)]] 
+ \Yj'^r[<PWK f i/r"5(r,r")V*,(^(r")V*,y(r",r)j 
+ zZ^^f^W r ^/r"V^,0(r")Vf„[5(r,r")V«y(r",r)]l 
-lYj^^'im [ flfr"5(r,r")Vf„^(r")Vf„V«y(r",r)] 

a b 



,D(4A2)(j.-) 

= ^^^Z/ ^(r'r")0(r')V^0(r")Vf„5(r,r") 

X 5(r, r') y(r", r') + ^(''^^'(r) 

= ^A|^Jfl!r"_^Vf„0(r")Vf„[,5(r,r")r(r",r)] 

- fl'r" ^ Vf„«^(r") Vf„r(r", r) 5(r, r" 



+ a)(4*4)(-r-) 



4 L-^ 2 
Combining together we obtain: 



= - ^ J] V^[vf 0(r)V° r Jr>(r')^(r, r') Vf y(r, r') 

aft ^ 

aft ^ 

a b 

$(463)(r) = - J] Vf^VCV^ • Vn) - 20 J] r^VVj, 

a b 

(D(4''4)(r) = - ^ r i/r"Vf„0(r")r(r",r)Vf„V^^(r,r")l 

aft ^ 

= -\Y,^'f[mK r rfr"V';„0(r")r(r",r) Vf„5(r,r")] 

aft ^ 

+ zZ^"kW r ^'r"V*^,0(r")V«r(r",r)Vf„5(r,r")l. 



<D(4fc)(r) = -iA[0(V0 • V«)] + ■ [0V(Vn • V0)] 

- ^Af^nA^l + iv • [^*V(nA<*)l 

- ^ V''[<^T''V*0] - ■ [(/.A0Vn]. (A17) 

a 

Equation for O'"''^' is obtained from (A 16) with replacement 
^(r,r") ^ r(r,r") and r(r",r') ^ 5(r",r'): 

0(4c)(r) = 2]^" r d{Y'r")<p{Y')V^^,<p{r") 

ab ^ 

X [y(r, r") V^5(r, r') Vf„5(r", r') 

- Vf„y(r, r") V^(5(r, r') ^(r", r') 

- V^7(r, r") 5(r, r') Vf„5(r", r') 

+ V*, V°(7(r, r")) S{v, r') J(r", r')]. (A18) 
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Integrations yield: 

ab 

3,(4c3)(r) = _1 ^ v°<^ r i;r"Vf„0(r") VXr,r")V*,J(r",r) 

It is little bit more difficult to get the first term: 
a,(4^i)(r) = - ^ a[,^ J dr" r(r, r") ^ V^,0(r") V^,5(r", r)] 
+ zZ^^^ r «?r"<-^(r")V«r(r,r")Vf„5(r",r) 

ab ^ 

= ^AUnA^l + ^aU(V0 • Vn)l 

4 L J g L J 

a ab 

Combining together we obtain: 

a)(4c)(-r) = 1a[0(V0 • Vn)] + ^A[0nA0] 

- ^ V''[(^T''*VV] - ■ [(^A0Vn]. (A19) 

a 

Finally substituting eqs. (A14), (A17) and (A19) in (A5) we 
obtain for the forth order term: 

^C^Xr) = i V • [V(^ (nA0 + (V0 • Vn)) + ^nV(V(^)^]. (A20) 

<l)(r) has the meaning of an induced density, i.e. it is a gradi- 
ent of the non-linear polarization. Expression for this quantity 
was derived by Wang, Cheng and Bower [26] and Apell [27] 
from the classical equation of motion for the electron in elec- 
tromagnetic fields. Our expression is identical to their result 
(cf. Eq. (11) in [26, 27]) if we write the external field in the 
form E = -V0. Our expression is more general in a sense that 
E also contains the induced fields. Note there are no more 
terms containing y{r,r") in the expression (A20). The high- 
frequency condition is not the only way to simplify Eq. (Al). 
In particular, Rudnick and Stem [39] demonstrated another 
possibility to arrive at the classical expression (Eq. (1)) start- 
ing from the quantum mechanical response function for the 
homogeneous electron gas and applying the local approxima- 
tion qvp (jj, where vp is the Fermi velocity and q is the 
inverse of the characteristic length of variation of the field. 



Appendix B: Non-linear source term ^^\r) 

Here we present expressions for the second part of the 
second-harmonic source term. Calculations are faciUtated by 
the use of maple computer algebra package. 



(2b) 



1 



6/1/2 



(m'-8n)+^(/i/2'+ 2/2/1') 



^(7/in2 + 15ni/2) - ^(5/1/2' + 9/2/;) 
+ ^^rn' - 4n)(ni/2' + "2/1' + ^) 



+ Ann[f[n'2 + f2n[ - 87rnin2) 



(Bl) 



47m{fin\-47ml) + ^(Sr/im -4/;ni) 



4w4V57r 

.n'(4.n,/(.M)!_^_|) 



(B2) 



We use the same abbreviated notations as in Sec. in and ad- 
ditionally introduce ^^''^ = n = n{r). 



Appendix C: Fast numerical solution of the integral equation 

It is convenient to recast eqs. (24) and (25) in the following 
form: 

.(0)/ 



G'r\r,r')ae{r';co), (CI) 



where for the linear response we have ^e(r, co) = [af\r, a)). 
Let us now omit the a> argument for simpification and use the 
definition (23): 

+ 1 



ae{r) = ^e{r) + af{r)[(a- ^«^i(r) - {w2{r)\, 



with 



!«iW = ^ {r'f\t{r')dr', 

W2(r) - j I — I a({r)dr , and a-W2{oo). 
Functions wi,2(?") satisfy following differential equations: 
w\(r) = /"-^^^({r) + (0 [(a - ^«'i(0 - (imir)^, 

"^2('-) = -^^\^t{r) + af{r) [{a - ^"^i(r) - ^«;2(r)]|. 

They can be solved at the 0(N) cost where A'^ is the number of 
mesh points to represent the density. 
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